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INTRODUCTION 


Fire  Control  system  needs  for  increased  probabilities 
of  detection,  recognition,  and  identification  require  the 
improvement  of  tactical  target  edge  and  feature  contrast, 
and  signal-to-noise  ratio,  i.e.,  image  enhancement.  Since 
optical  techniques  can  provide  only  limited  algorithmic 
flexibility  and  image  enhancement,  more  powerful  digital 
techniques  such  as  spectral  filtering  must  be  tested. 

The  work,  begun  in  1975,  relating  to  the  development 
of  a  digital  image  enhancement  interactive  demonstration 
program  is  presented  in  this  paper.  Originally  existing 

on  the  West  Point  U.S.  Military  Academic  computer,  it  now 
exists  on  the  U.S.  Army  ARRADCOM  CDC  6500  time-sharing 
system. 

The  motivation  for  this  effort  lies  in  work  done  2  years 
previously  by  a  team  of  investigators  at  Northrup  Research 
Corporation  in  Anaheim,  California.  The  researchers  did  obtain 
some  image  enhancement  by  using  a  hard-wire  Haar  transform 
processor  for  contrast  improvement.  Since  these  results 
were  inconclusive,  the  decision  was  made  to  carry  out  an 
in-house,  laboratory,  independent  research  effort  with  the 
objective  of  using  digital  image  enhancement  techniques 
for  contrast  and  signal-to-noise  improvement  as  an  aid  in 
target  detection,  tracking,  and  weapon  direction.  The 
currently  existing  routine,  Program  IMAGE,  is  an  outgrowth 
of  this  effort. 

First,  some  mathematical  background  is  presented  in 
matrix  theory,  Fourier  analysis,  and  Hadamard  matrices. 

Next,  the  steps  taken  in  the  development  of  Program  IMAGE 
and  an  explanation  of  the  roles  of  and  tasks  performed  by 
the  various  modular  subroutines  are  discussed.  Finally, 
several  examples  of  actual  image  enhancement  computer 
outputs  are  presented  and  explained . 

MATHEMATICAL  BACKGROUND 
Matrix  Manipulation 

Arbitrary  Linear  Gray-Scale  Transformations 

Gray-level  scaling,  the  capacity  to  reset  the 
digitized  values  of  the  individual  picture  elements  of  an 
image  to  within  new  upper  and  lower  bounds,  is  important 
for  two  reasons.  First,  after  algorithmic  processing,  it 
may  be  important  to  rescale  the  output  data  to  put  the 
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bounds  within  meaningful  limits  for  easy  interpretation, 
spotting  of  anomalous  results ,  and  the  generation  of 
histograms  that  can  be  quickly  and  easily  compared  with 
those  of  the  original  input  data.  Second,  it  is  important 
to  place  the  lower  bound  of  image  processing  results  just 
above  the  minimum  level  which  is  available  on  a  digital  to 
analog  video  display  screen,  to  prevent  "clipping"  of  the 
minimum  values.  Similarly,  the  upper  data  bound  should 
be  placed  below  the  maximum  gray  level  setting.  This 
prevents  the  "whiting"  or  "washing  out"  of  the  results,  or 
even  worse,  the  loss  of  valuable  picture  information. 

Problem:  Given  that  picture  element  values  range 

between  A  and  B,  rescale  them  to  range  between  the  values 
C  and  D.  That  is,  for  A  S  X  g  B,  find  the  linear 
transformation  function,  y(x),  such  that  C  §  y  s  D. 

and : 


y  (A)  =  C 

y  (B)  =  D  {1) 

From  the  familiar  point-slope  equation  for  a  linear 
function,  for  a  given  point  (x-,,  y-,)  and  any  point  (x,  y) 
one  has : 

y  - .  yx  =  M  (x  -  X-.) 

1  1  (2) 

M  =  D  -  C 

B  -  A  (3) 

Take  xp  =  A,  and  yp  =  C.  If  one  took  x-!  =  B,  and 

yp  =  D,  the  same  end  result  would  be  achieved.  Substituting 
for  Xp  and  yp  into  Equation  2 , 

V  -  C  =  D  -  C  •  (x-A) 

B~-~ A  (4) 

Adding  C  to  both  sides  under  the  denominator  B  -  A  and 
rearranging  terms  slightly  results  in  the  arbitrary  linear 
gray-scale  transformation  equation, 

y  =  D_ C  .  x  -  A'D  -  BC  .  . 

B  -  A  B  -  A 


By  this  equation,  pixel  data,  ranging  from  values  A  to  B, 
are  linearly  transformed  to  a  set  of  values  ranging  from 
C  to  D. 


As  a  simple  example,  one  can  calculate  the 
transformation  scaling  gray-level  values  ranging  from 
0  to  1  into  gray-level  values  ranging  from  0  to  100. 
In  this  case,  substitute  into  Equation  5: 

A  =  0 


B  =  1 
C  =  0 


(6) 


D  =  100 

As  one  might  expect,  in  this  instance,  Equation  5  simply 
reduces  to: 

y  =  lOOx.  (6) 


Original  Data  Scan  and  Commutativity 


As  mentioned  in  the  introduction,  motivation  for 
this  work  was  derived  from  the  Northrup  Corporation , 
specifically,  from  a  one-dimensional  Haar  transform  scan 
that  they  used.  The  scan  forms  a  transform  image  whose 
rows  are  formed  by  the  operation  of  the  rows  of  the  Haar 
matrix  T,  on  the  rows  of  the  image  matrix  I.  This  is  in 
contrast  to  standard  matrix  multiplication  where  the  rows 
of  the  product  matrix  are  formed  by  operating  with  the  rows 
of  the  transformed  matrix  on  the  columns  of  the  matrix  of 
the  image  to  be  transformed  and/or  enhanced.  With  a  simplified 
one-subscript  notation  for  easier  reading,  one  can,  for 
example,  multiply  the  2-by-2  Haar  transform  matrix  T,  by  the 
rows  of  the  image  matrix  I.  To  obtain  the  function  F(T,  I), 
a  formal  expression  for  the  representation  of  the  image 
transform,  one  has: 
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F(T,  I)  =  T  -IT  (8) 

T 

where  I  is  the  transpose  of  the  image  matrix  I.  This  is 
required  so  that  the  rows  of  the  transform  matrix  T 
multiply  the  rows  of  the  image  matrix  I,  and  not  the 
columns.  The  result,  for  the  one-dimensional  horizontal 
scan  is : 


F  (T,  I)  = 


T1  xi  +  t2  x2  t1  13  +  t2  *4 


t3  ix  +  t4  i2  t3  i3  +  i4  I, 


(9) 


In  matrix  manipulation,  the  question  arises 
as  to  when  the  product  of  two  matrices  is  commutative, 
i.e.,  when  A  B  =  B  A,  for  two  N-by-N  matrices  A  and  B. 
For  the  2-by-2  case: 


(10) 


Take  the  first  component  of  the  two  product  matrices,  for 
example : 


or 


A1  +  A  2  B3 

a2  =  A3  ^2 


(ID 


One  sees,  therefore,  by  comparing  with  other  components 
that  A  •  B  =  B  •  A  when  the  matrices  are  symmetric. 

Two-dimensional  Scanning  and  the  Role  of  the  Multipliers 

Given  an  image  matrix  I,  a  transform  matrix  T,  an 
inverse  transform  matrix  T--*-,  an  output  matrix  O,  and  a 
storage  matrix  X,  one  has  the  relationships  (ref.  1): 

X  =  T  •  I  •  Tt 

(12) 

0  =  T_1  •  X  •  (T“1)T 
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The  Haar  matrix  T,  in  the  4-by-4  case,  contains 
values  or  enhancement  factors  (ref.  2)  .  It  is  of  the  form 


T  = 

M1 

Mi 

Mi 

M1 

m2 

m2 

M2 

M2 

m3 

m3 

M3 

m3 

m3 

m3 

m3 

M3 

(13 

The 

.  T 

matrix  T  , 

is  proportional  to 

a  matrix  whose  columns 

are 

identical 

to  the  rows 

of  the  one  in  Equation 

13. 

The 

matrix  of 

the 

image  is  of 

the  form: 

I  = 

r 

xii 

x12 

I13 

Tu 

*21 

I22 

I23 

I24 

*31 

I32 

I33 

Z34 

Z41 

Z42 

Z43 

I44 

(14: 

To  see  how  the  M's  behave  under  transformation, 
that  the  storage  matrix  X  is  proportional  to  T*I-Tt. 


Performing  the 

product  I 

•  tt 

yields 

for : 

si  =  in 

+ 

Z12 

+ 

I13 

+ 

I14 

s2  =  121 

+ 

122 

+ 

I23 

+ 

I24 

(15a) 

S3  "  I31 

+ 

I32 

+ 

x33 

+ 

X34 

S4  =  J41 

+ 

J42 

+ 

143 

+ 

J44 

S1M1 

S1M2 

S1M3 

S1M3 

S2Mi 

®2^2 

S2m3 

s2M3 

S3M1 

S3M2 

S3M3 

S3M3 

S4Mi 

S4M2 

S4M3 

S4M3 

three  M 


note 
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Continuing  the  procedure  and  multiplying  the  matrix  expression 
in  equation  15b  by  the  matrix  T  on  the  left  result';  in 


X  =*=  F  »  M  (16) 

where  F  depends  on  the  sum  of  the  image  data  values, 

F=I11+i12+i13+I14+I21+i22+i23+i24+i31+i32+i33+i34+*41+i42+I43+i44 

(17) 


and  M  is  a  matrix 

,2 


M  = 


m: 


Mi  M2 


M1  M3 


Mi  M3 


Mx  M2 


m; 


m2  m3 


m2  m3 


M!  M3 


m2  m3 


Mi  M- 


m2  m3 


m: 


m 


3 

2 


M_ 


(18) 


Thus,  two-dimensional  image  transform  matrices  are 
quadratically  dependent  upon  the  enhancement  factors  in  a 
symmetric  matrix  format,  with  the  main  diagonal  being 
quadratic  in  terms  of  purely  one  kind. 

Fourier  Transform  Background 

Summation  and  Matrix  Equivalency 


A  given  N-by-N  digital  image  matrix  has  picture 
element  entries  f  where  the  indices  x  and  y  indicate  the 
row  and  column  positions  and  are  taken  to  range  from 
0  to  N-l.  Its  Fourier  transform  is  the  N-by-N  matrix 
with  elements  Fuv,  and  is  given  by  the  familiar  two- 
dimensional  Fourier  transform  expression: 


F 

uv 


N-l  N-l 


/'-2  m  i 
exp  \  N 


(ux  + 


(19) 


For  the  matrix: 


1 

1 

1 

1 

1 

w-1 

-2 

w  •  •  • 

w-(N-l) 

1 

w-2 

w-4  ... 

w-2(N-1) 

. 

. 

• 

• 

• 

• 

• 

• 

• 

• 

1 

-2 (N-l) 
w  ... 

W-<N-1)2 

where 


(20) 


w  =  e 


.2  jt  i 
N 


(21) 


One  can  write  a  matrix  equation: 

F  =  W  •  f  •  W  (22) 

Problem:  Show  that  Equations  19  and  22  are 

equivalent.  Take,  as  an  example,  N  =  3,  and  calculate 
the  first  element  Fqq.  From  Equation  19,  it  is  evident  that 

F00  =  f00+f01+f02+fl0+fH+f12+f20+f21+f22  (23) 


In  this  case,  from  the  definition  of  f,  and  Equation  20 


To  compute  Fqq  by  means  of  Equation  22,  take  the  matrix 
X  =  F * W  and  compute  the  first  column  of  X 


xll~ 

foo 

+ 

f01 

+ 

f  02 

X21= 

f10 

+ 

fll 

+ 

f  12 

(25) 

X31 

=f  20 

+ 

f  2 1 

+ 

f  22 

To  obtain  Fqq,  multiply  X  on  the  left  by  the  first  row  of 
the  Fourier  transform  matrix  of  W  of  Equation  20. 

F00  =  f00+f0i+f02  +  f10+fU+f12+f20+f21+f22  (26) 

This  expression  is  identical  to  Equation  23.  Thus,  the 
summation  and  matrix  notations  are  equivalent. 

Delta  Function  and  Normalization 

Image  enhancement  can  work  because,  from 
sampling  theory,  transforms  can  select  out  and  increase 
the  prominence  of  a  given  spatial  frequency  component  of 
a  signal.  The  Nth  spatial  frequency  component  of  a  signal 
is  defined  as  N  variations  from  a  low  to  a  high  pixel  value, 
over  the  data  field  length.  The  delta  function  has  the 
property  of  being  able  to  select  out  a  given  value. 
Therefore,  the  derivation  of  its  integral  expression  will 
be  reviewed  and  a  corresponding  summation  expression 
derived  for  application  to  digital  data  sets. 


For  a  function  of  time,  f(t),  and  its  Fourier 

transform  ,  F(w) 


/i  wt 
F  ( w)  e  d  w 

F (w)  =  J  f(T)e“lw7dT 


Substituting  Equation  28  into  Equation  27  and 
reversing  the  order  of  integration  because  the  functions 
are  assumed  to  be  continuous  and  integrable 

r  jw(t  -r)  l 


f (t)  =  /fcn[  j 

If  one  writes 

f  (t)  =  j  f  (T)  5  (  T-T)  d' 


(T-T)  d 7 


P 


In  this  case,  the  delta  function  filters  out  the  value 
f(t-T)  to  be  input  into  f{t).  Then  one  has 


5  (t-T  )  =  fz  ±w{t~T)dT 


(31) 


This  is  the  integral  expression 
Similarly,  for  discrete  digital 
digital  function  of  the  index  x 
equivalent 


for  the  delta  function. 

functions  where  fx  is  the 

and  Fjj  is  the  Fourier  transform 


f 


x 


where  w 


N-l  F  w  U  x 
(1=0  ^ 

e^  31  as  before.  Also, 


y~N_1  fxi  w_^?xi 
xi=0 


(32) 


(33) 


By  substituting  Equation  33  into  Equation  32  and 
reversing  the  order  of  the  summation  indices: 


And,  therefore,  as  in  Equation  31: 

us) 

The  last  point  to  consider,  with  regard  to  Fourier 
transforms,  is  the  question  of  normalization  for  proper 
scaling,  a  point  not  always  explicitly  treated  in  the 
literature.  As  an  example,  take  Equation  19,  for  the 
complex  (real  and  imaginary) ,  two-dimensional  Fourier 
transform  and  reduce  it  to  the  one-dimensional  case  for 
the  real  part  of  the  Fourier  transform  of  a  real  function: 


f  j  Cos 


(36) 


where,  from  sampling  theory,  the  permissible  values  of  i 
are,  in  general,  0  through  N-l. 

2 

Take,  as  an  example,  the  case  of  N  =  8 .  If  the 
function  of  the  Fourier  transform  that  is  being  sought  is 
one  cycle  of  a  cosine,  then  Fg,  the  DC  term,  is  0  and  F^, 
the  weight  of  the  fundamental  cosine  frequency  is  required 
to  be  equal  to  1.  Substituting  Fi  =  1  into  Equation  36  and 
dividing  by  a  constant  of  proportionality  N  for  scaling 
yields  s 


9 


1 


-zu  « 


(37) 


8 

s 


Here,  s  is  the  unknown  scaling  factor  and  gives 


f  j  -Cos 
Therefore , 


2  jt  j 
N 


(38) 


1  = 


i 


j=0  Cos' 


2  j 


(39) 


8 

s 


This  equation  yields: 
s  =  2. 


(40) 


Thus,  for  proper  scaling,  the  results  of  a  one-dimensional 
Fourier  transform  must  be  divided  by  N  and,  as  can  similarly 

2 

be  shown,  those  of  an  inverse  Fourier  transform  must  be 
divided  by  2.  They  must  share  a  common  divisor  of  N. 
Similarly,  the  results  of  a  two-dimensional  Fourier 
transform  must  be  divided  by  and  its  inverse  by  4. 

4  _ 


Therefore,  a  common  divisor  of  Nz  is  shared. 
Iladamard  Transform  Matrix 


The  well-known  Hadamard  matrices  are  based  on  the 
arrangement  of  Walsh  functions.  The  first  and  simplest  is 
the  2-by-2  case: 


H1  = 


1 

-1 


(41) 


:  o 


ri'i  i  i  ii 


Matrices  of  higher  order  are  recursively  generated  or 
"grown"  by  the  replacement  of  each  element  of  the  matrix 
in  Equation  41  by  the  matrix  H^ .  This  means,  of  course, 
that  when  the  -1  is  replaced  by  Hi ,  the  signs,  obviously, 
are  simply  reversed.  The  next  Hadamard  matrix,  the  4-by-4 
case  is 


Note  that  the  matrix  is  symmetric  and  except  for  a  factor 


of  4,  H  is  identical  to  H  for  any  order.  Here, 

=  1  .  H,  (43) 

2  4 

This  recursive  process  can  be  continued  indefinitely. 

The  N-by-N  matrix  is  multiplied  by  the  N  element  pixel 
vector,  and  the  elements  of  the  output  column  transform 
vector  are  adjusted.  For  example,  by  decreasing  the 
magnitude  of  the  first  element,  d.c.  or  background  energy 
is  suppressed.  The  inverse  transform  matrix  is  then 
applied,  and  the  result  is  the  output  pixel  vector. 

PROGRAM  IMAGE 

Development  of  the  Program 

Two  16-by-16  image  matrices  for  each  printout  page 
are  used,  the  upper  one  for  the  low-contrast  input  image, 
the  lower  one  for  the  improved  contrast  algorithm-processed 
output  image . 

The  elements  in  the  input  and  output  target  image  of 
the  simulated  image  matrices  correspond  to  the  brightness 
of  points  (or  areas  of  points)  of  a  larger  actual  TV  image. 
The  targets  are  characterized  by  a  bordering  background 
level  and  a  central  signal  level  outlined  by  a  graphical 
drawing  of  a  tank.  The  characterization  of  the  tank  target 
matrix  is  shown  in  figure  1.  The  letters  G,  C,  T,  H,  W,  and 
B  show  the  various  sections  of  the  tank  in  the  image  field 
and  indicate  the  gun,  cupola,  turret,  hull,  the  wheels 
(or  tread),  and  the  background,  respectively.  (Appreciation 
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Figure  1.  Tank  target  character! zation, 


is  extended  to  a  former  co-worker  Gerald  VTatson  for  his 
assistance  in  implementing  the  computer  programming  and 
the  computer  graphics  required  for  a  meaningful  display 
of  the  output  data.) 

No  matrix  larger  than  the  16-by-16  case  has  been  used 
in  this  effort  for  two  reasons:  First,  the  fast  Haar  and 
the  fast  Fourier  transforms  used  algorithms  that  require 
the  dimension  N  of  the  input  data  to  be  an  integer  power 
of  two.  Thus,  the  next  larger  utilizable  square  image 
matrix,  the  32-by-32  case,  contains  1024  entries  and  is 
too  unwieldy  to  easily  manipulate  on  paper,  both  for 
theoretical  studies  and  numerical  value  checking.  Also, 
the  16-by-16  image  matrix  is  the  largest  one  that  can  be 
fully  displayed  on  the  computer  terminal  because  of  con¬ 
straints  in  the  number  of  alphanumeric  characters  and 
spaces  available  per  output  line. 

The  chronological  development  of  Program  IMAGE  is 
summarized  in  table  1.  The  first  option  programmed  from 
equations  supplied  by  the  Northrup  investigators,  w as  the 
Haar  4-by-4  matrix  that  affords  the  capacity  to  vary  three 
spatial  frequencies  for  image  enhancement,  the  highest  being 
value  changes  between  adjacent  picture  elements.  After 
initial  tests,  another  processor,  a  Haar  16-by-16  matrix, 
was  added  that  affords  greater  flexibility  via  the  capacity 
to  manipulate  two  more  spatial  frequencies  for  image 
enhancement.  Specifically,  a  d.c.  term,  and  1,  2,  4,  and 
8  cycles  of  square  wave  oscillation  are  available  for  the 
data  field  length  of  16  points.  For  a  more  detailed 
discussion  of  this,  see  figure  3  and  the  relevant  section 
on  matrices  in  the  referenced  Haar  transform  report  (ref.  2) 

Only  a  background  and  signal  level  were  initially 
displayed  to  get  the  software  working  as  quickly  as  possible 
Later  a  random  noise-level  option  was  added  to  give  a  more 
real-world  situation  whereby  one  could  process  TV  pictures 
with  moderate  or  severe  noise  levels  incorporated  into  the 
data . 


To  chart  how  signal-to-noise  ratios  propagate  through 
algorithms  and  to  see  if  they  can  be  improved,  a  subroutine 
was  added  to  the  program  to  calculate  the  RMS  value  of  the 
processed  picture  patterns.  Since  transforms  or  other 
algorithms  enhance  or  pull  out  contrast  either  horizontally 
across  a  line  thus  causing  vertical  edges  to  stand  out,  or 
down  a  vertical  column  making  horizontal  lines  stand  out, 
or  in  two  dimensions,  an  option  was  added  to  calculate  and 
print  out  contrasts  for  each  section  of  the  image  in  both 
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Table  1.  Development  of  program  IMAGE 


Steps 

1 


3 

4 

5 

6 

7 

8 

9 

10 
11 


Item 

Establish  computer  formatting  and  programming. 

4-by-4  Haar  transform  matrix,  (Northrup  Corp's 
suggestion) . 

16-by-16  Haar  transform  matrix. 

Random  noise  capability. 

Processed  RMS  value  of  image,  capability. 
Contrast  calculations. 

Average  processed  signal  option. 

Modularization . 

Subroutine  BEEP. 

Corrections  to  RMS  noise-level  scaling. 

Current  status  of  nine  available  algorithms. 


horizontal  and  vertical  directions.  Appreciation  is 
extended  to  ARRADCOM  engineer,  Mr.  Ray  Popko,  for  his 
helpful  suggestion  of  this  idea. 


A  routine  employing  the  4-by-4  Hadamard  transform 
in  a  horizontal  scan  was  written  and  separately  pro¬ 
grammed.  For  easier  access,  it  was  added  to  the  already 
existing  Haar  4-by-4  matrix  and  the  Haar  16-by-16  matrix 
routines,  and  combined  into  one  file.  In  an  effort  to 
provide  even  more  flexibility  in  characterizing  and 
manipulating  the  spatial  frequency  components  of  image 
data  for  enhancement  purposes,  a  Fourier  transform  test 
program  was  written  to  understand  its  numerical  behavior 
in  one-  and  two-  dimensional  scans. 

The  Fourier  transform  provides,  for  N  points,  N 

2 

spatial  frequency  values  that  can  be  weighted  for  enhance¬ 
ment.  Upon  debugging,  Fourier  transform  and  Fast  Fourier 
transform  routines  were  incorporated  in  horizontal  and 
two-dimensional  scans. 


'  £ 


An  option  was  then  added  to  compute  the  average  processed 
signal  strength  of  the  tank  target  in  the  data  field  for 
signal-to-noise  measurement  purposes. 

Because  of  the  increasing  difficulty  in  modifying 
the  programming/  then  stored  on  the  West  Point  Honeywell 
computer,  the  decision  was  made  to  rewrite  the  entire 
program  in  terms  of  transform-computing  and  picture¬ 
processing  modular  subroutine  packages  that  could  simply 
be  strung  together.  At  this  time,  the  cumbersome  matrices 
and  equations  of  the  Haar  4-by-4  and  16-by-16  routines  were 
replaced  by  the  author's  fast-Haar  transform  algorithm 
that  requires  no  two-dimensional  matrix  arrays .  This 
action  reduced  the  bulk  of  the  coding  and  reduced  the 
storage  area  requirements  of  the  program  by  60  percent 
(ref.  2) . 

Lastly,  a  more  detailed  signal-to-noise  analysis  was 
conducted,  resulting  in  the  proper  scaling  of  noise  values 
such  that  their  RMS  value  becomes  equal  to  a  desired  value 
without  error.  This  derivation  is  discussed  in  full  detail 
in  the  section  on  RMS  noise-level  scaling  in  the  referenced 
report  on  the  application  of  the  author ’ s  fast-Haar  transform 
algorithm  to  real  IR  data  (ref.  3). 

Utilizing  the  Program 

The  divisions  of  the  total  main  program/subroutine 
package  and  IMAGE  program  are  outlined  in  table  2.  The 
modular  subroutine  packages  are  arranged  in  the  order  in 
which  they  are  called  by  the  main  program  with  Haar 
followed  by  Hadamard,  followed  by  Fourier  transform 
processing  routines. 

A  summary  of  the  current  status  of  the  program 
involving  nine  presently  available  algorithm  options  for 
scanning  imagery  data  in  horizontal,  vertical,  and  two- 
dimensional  formats  is  given  in  table  3. 

The  numerical  data  in  the  beginning  of  the  program 
(see  appendix  A)  define  the  coordinates  of  the  lines 
comprising  the  outline  of  the  tank  target  area.  After 
this  part,  the  program  requires  the  answers  to  one  question 
and  two  commands . 

1.  Do  you  want  to  test  number  sequences  on  your  algorithm? 

2.  Enter  background,  signal,  and  noise  level. 

3.  Specify  which  of  the  following  algorithms  by  number. 
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Table  2.  Division 
No.  Section 

1  Main  program 

2  Subroutine  BEEP 

3  Subroutine  APSIG 

4  Subroutine  PRMS 

5  Subroutine  CTRAST 

6  Subroutine  HAAR 

7  Subroutine  HAR16H 

8  Subroutine  HAR16V 

9  Subroutine  HAR16HV 

10  Subroutine  HAR4H 

11  Subroutine  HAD 

12  Subroutine  HAD4H 

13  Subroutine  FT 


of  program  IMAGE 
Function 

Signal  parameter  input, 
alphanumeric  output, 
graphics,  calls  to  subroutines. 

Prompts  user  when  ready  for 
input . 

Calculates  average  processed 
signal  level. 

Calculates  RMS  value  of 
processed  image. 

Edge  contrast  of  features  of 
tank  target. 

Investigator's  fast  Haar 
transform  algorithm. 

Effects  horizontal  scan 
using  subroutine  HAAR. 

Effects  vertical  scan  using 
subroutine  HAAR. 

Effects  horizontal  then 
vertical  scan,  using 
Subroutine  HAAR. 

Horizontal  scan,  in  4-pixel 
blocks,  using  Subroutine 
HAAR. 

4-pixel  Hadamard  transform. 

Horizontal  scan  in  4-pixel 
blocks,  using  Subroutine 
HAD. 

Effects  conventional,  one¬ 
dimensional  Fourier 
transform. 


14 


Subroutine  FT16H 


Horizontal  scan  using 
Subroutine  FT. 


Table 

2.  Division  of  program 

IMAGE  (continued) 

No. 

Section 

Function 

15 

Subroutine  FT2D 

Conventional  two-dimensional 
Fourier  transform. 

16 

Subroutine  FT16D2 

Effects  two-dimensional 
Fourier  transform  of 

16-by-16  image. 

17 

Subroutine  FFT  and 
REVERS 

Implements  one-dimensional 
Fast  Fourier  transform. 

18 

Subroutine  FFTH 

Effects  horizontal  scan 
processing  by  the  Fast 
Fourier  Transform/algorithm. 

19 

Subroutine  FFTHV 

Effects  horizontal  then 
vertical  scan  processing 
using  the  Fast  Fourier 
Transform  algorithm. 

Table  3.  IMAGE  processing  options 

No, 

Transform  type 

Processing  procedure 

1 

Haar  16-  element 

Horizontal  scan 

2 

Haar  16-  element 

Vertical  scan 

3 

Haar  16-  element 

Horizontal  then 
vertical  scan 

4 

Haar  4-  element 

Horizontal  scan 

5 

Hadamard  4-by-4 

Horizontal  scan 

6 

16-  element  Fourier 

Horizontal  scan 

7 

16-by-16  grid  Fourier 

Two-dimensional  scan 

8 

16-  element  Fast  Fourier  Horizontal  scan 

9 

16-  element  Fast  Fourier  Horizontal  then 

vertical  scan 
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If  one  replies  "yes"  to  the  "number  sequence"  option, 
the  16-by-16  input  image  matrix  array  is  bypassed.  One 
can  then  enter  either  a  one-dimensional  test  signal  vector 
for  algorithms  one  and  five,  a  16-element  input  to  the  Ilaar 
transform  algorithm  or  a  4-element  input  to  the  Hadamard 
transform  routine,  respectively.  The  resulting  16-  or  4- 
element  series  is  printed  out  after  processing  by  enhance¬ 
ment  weights  that  the  user  selects.  If  the  number  sequence 
option  is  not  invoked,  the  user  is  required  to  furnish  the 
background  signal-and-noise  level  and  a  given  transform. 
Suggested  data  selections  that  one  can  use  to  test  for 
improved  contrast  with  initial  input  contrast  ranging  from 
5  to  50  percent,  every  5  percent,  are  listed  in  table  4. 

Similarly,  if  one  is  primarily  interested  in  signal- 
to-noise  analysis,  suggested  input  parameters  for  input 
signal-to-noise  ratios  ranging  from  1.5  to  6.0  every  0.5 
unit  are  listed  in  table  5. 

Upon  assignment  of  the  desired  algorithm,  the  program 
will  indicate,  if  needed,  the  enhancement  weighted  factors 
to  be  tested.  As  will  be  shown  in  the  examples  later, 
good  enhancement  will  be  achieved  by  the  first  algorithm, 
a  16-element  Haar  transform  scan  with,  for  example,  the 
values  0.5,  ].0,  5.0,  10.0,  and  15.0. 
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Table 

4. 

Suggested  data  selections 

for  contrast  demonstrations 

Background 

Signal 

Noise 

Signal  contrast 

(%) 

19 

20 

0 

5 

9 

10 

0 

10 

17 

20 

0 

15 

8 

10 

0 

20 

9 

12 

0 

25 

14 

20 

0 

30 

13 

20 

0 

35 

12 

20 

0 

40 

11 

20 

0 

45 

10 

20 

0 

50 

Table 

5. 

Suggested  input 

parameter 

selections  for  signal- 

. 

to-noise  testing 

Background 

Signal 

Noise 

Signal-to-noise 

ratio 

5 

9 

6 

1.5 

5 

10 

5 

2.0 

6 

10 

4 

2.5 

5 

9 

3 

3.0 

8 

14 

4 

3.5 

8 

16 

4 

4.0 

10 

18 

4 

4.5 

6 

10 

2 

5.0 

6 

11 

2 

5.5 

6 

12 

2 

6.0 
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Tasks  of  the  Modular  Subroutines 

Now,  a  brief  discussion  follows  in  order  of  the 
coding  of  the  various  modular  subroutines  and  the  tasks 
by  which  they  are  performed.  Each  algorithm  is  called 
from  the  main  program  by  the  parameter  "INDEX"  whose 
legal  values  are  integers  ranging  from  one  through  nine. 
With  the  addition  by  the  user  of  more  image  processing 
algorithms  to  the  program,  this  can  easily  be  modified. 

Subroutine  BEEP 

This  routine  prompts  the  user  with  a  series  of 
audible  bells  as  to  when  the  terminal  is  ready  for  data 
input.  The  three  call  parameters  represent,  in  turn, 
the  number  of  bells,  the  length  of  time  (approximately  in 
seconds)  of  each  bell,  and  the  length  of  the  delay  time 
(approximately  in  seconds)  between  each  bell.  The  best 
results  assume  a  CDC  6500  system  tied  through  a  SCOPE 
operating  system  to  a  Tektronix  4014  graphics  terminal. 
Bell  sounding  and  blank  delay  time  characters  may  pass  at 
different  rates  on  other  terminals.  Also,  on  other 
computing  machines  and/or  operating  systems,  the  octal 
definitions  of  the  bell  and  blank  variables  may  not 
necessarily  apply. 

Subroutine  APSIG 

This  subroutine  calculates  and  prints  out  the 
average  value  of  the  strength  of  the  processed  signal. 

The  pixel  values  for  the  tank  target  are  simply  added  for 
the  gun  and  cupola,  turrent,  hull,  and  wheels  in  four 
do  loops.  The  average  is  then  taken  and  the  resulting 
value  printed  out. 

Subroutine  PRMS 

This  subroutine  calculates  the  RMS  value  of  the 
processed  output  image  and  is  useful  for  testing  noise 
level  propagation  under  the  influence  of  algorithms  of 
interest.  All  of  the  pixel  values  over  the  entire 
processed  output  image  field  are  squared  and  added, 
averaged,  and  then  the  square  root  is  taken  before  print¬ 
out  . 

Subroutine  CTRAST 

This  subroutine  calculates  the  signal  contrast 
of  various  parts  of  the  tank,  i.e.,  gun,  cupola,  turret, 


hull,  tread,  and  finally  an  overall  average  value  for  the 
entire  tank  target  image  in  both  vertical  and  horizontal 
directions.  The  definition  chosen  for  output  signal 
contrast  is: 


C  =  S-b 

S  (44) 

Where  S  is  the  signal  level  of  a  tank  pixel  and  B  is  the 
background  value  of  a  neighboring  element  not  on  the  tank 
but  in  the  surrounding  background/clutter  region.  This  is 
done  for  all  occurrences  of  signal  and  background  adjacent 
element  pairs  for  each  element  of  each  part  of  the  tank  in 
both  horizontal  and  vertical  directions. 

For  example,  take  S  to  be  the  value  of  the  13th 
element  from  the  left  in  the  third  row  down,  element  (3,  13), 
one  of  the  gun  elements.  Take  b1  to  be  the  value  of  background 
element  (3,  12)  and  b2  to  be  the  value  of  background  element 
(3,  14) ,  the  elements  to  the  left  and  right  of  the  given  gun 
element.  The  horizontal  contrast  at  the  (3,  13)  element,  the 
first  element  on  the  gun  tip  is  taken  to  be: 

S-bi  S-b2 

c  =  ~  *  ~  -  (45) 

2 

for  the  horizontal  contrast  of  the  first  gun  element.  This 
procedure  is  implemented  for  each  section  of  the  tank  in 
both  horizontal  and  vertical  directions.  The  values  are 
expressed  as  percentages  for  printout.  The  Hollerith 
headers  GUNH ,  GUNV,  TURH ,  TURV,  etc.,  indicate  the  type  of 
feature  and  contrast  specified,  for  example,  GUNH  indicates 
gun  horizontal  contrast  value. 

Subroutine  HAAR 

This  subroutine  implements  the  author's  fast 
Haar  transform  algorithm  (ref.  2).  The  input  string  of 
pixels  is  entered  by  the  user  into  the  first  N  of  2N  -  2 
positions  of  the  Q  array,  the  log  to  the  base  2  of  N  is 
supplied  by  the  user  to  the  variable  NU,  and  the  NU  +  one 
enhancement  factors  are  read  by  the  user  into  the  E  array. 

The  N  output  enhanced-pixels  are  returned  by  the  routine 
to  the  last  N  of  the  2N  -  2  positions  of  the  Q  array. 
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Haar  Transform  Subroutines 


HAR16H ,  HAR16V,  HAR16HV,  and  HAR4H  subroutines 
employ  subroutine  HAAR  to  process  the  16-by-16  input  image. 
The  first  scans  16  elements  at  a  time  horizontally;  the 
second,  vertically;  and  the  third,  horizontally  then 
vertically.  The  fourth  breaks  up  each  16-element  horizontal 
picture  line  into  4-element  sections  and  invokes  subroutine 
HAAR  to  process  the  lines,  4  elements  at  a  time. 

Subroutine  HAD 

This  subroutine  invokes  processing  by  a  4-by-4 
Hadamard  matrix  as  defined  earlier.  Although  the  call 
parameters  are  expressed  in  general  terms,  in  this  case, 

Nu  must  be  two  so  that  2nU  that  equals  N  is  equal  to  four. 
Four  enhancement  weighted  factors  E  are  specified,  and  four 
input  picture  elements  to  be  processed  are  sent  into  the 
Q  array. 

Subroutine  HAD4H 


This  subroutine  invokes  subroutine  HAD  to  do 
Hadamard  transform  processing  on  the  input  image.  Each 
horizontal  input  picture  line  is  broken  into  four  equal 
length  sections  and  processing  is  done  four  elements  at  a 
time  on  the  input  image . 

Subroutine  FT 


Where  N  is  2nu  for  positive  integral  values  of 
NU,  this  subroutine  takes  the  Fourier  transform  of  N 
elements  and  incorporates  any  alterations  by  n  enhancement 

2 

factors  present  in  the  E  array.  The  inverse  Fourier  transform 
is  then  taken,  and  the  new  pixel  values  returned  to  the 
real  input  array  XR.  Assume,  for  example,  N  =  8.  Any 
given  Fourier  component  F^  is  of  the  form: 


F  = 
o 


E 


f . 
3 


tr  ig 


(46) 


j=0 

where  trig  is  a  trigonometric  function  (sine  or 
cosine)  of  the  given  argument.  Ordinarily,  to  take  the 
Fourier  transform  of  N  data  elements,  time-consuming 
trigonometric  computations  are  required. 


When,  for  example,  i  =  3  and  j  =  7,  trigono¬ 
metric  functions  of  the  argument  21 n  are  required.  However, 

8 

the  same  answer  is  obtained  because  of  periodicity  by 
subtracting  2N  or  16  from  the  argument  and  using  5  n 

8 

instead.  Thus,  the  trigonometric  values  of  the  arguments 
2  k  *  1  to  2  Jt  •  (N-l)  can  be  computed  once  and  stored  in  a 
N  N 

"lookup"  table.  By  simply  subtracting  integral  multiples 
of  N  from  the  trigonometric  arguments  in  the  ter.us  of  the 
Fourier  transform,  the  stored  results  in  the  lookup  table 
can  be  used  to  select  out  the  proper  trigonometric  values 
for  any  combination  of  the  indices  i  and  j . 

Also,  for  N  =  8,  the  Whitaker  sampling  theorem, 
allows  four  frequency  components  of  a  signal  to  be  represented 
by  this  Fourier  transform,  specifically  a  dc  term  and  1, 

2,  and  3  cycles  of  sinusoidal  oscillation  for  the  data 
field  length  of  8  points.  Just  four  frequencies  are  allowed 
because  at  least  two  samples  are  required  to  fully  represent 
a  signal  frequency  component.  However,  since  the  computed 
Fourier  transform  array  xr  contains  8  values,  it  was  found 
by  experimentation  that: 

xr^  =  dc  term 

xr2  =  strength  of  first  frequency 
xr^  =  strength  of  second  frequency 

xr4  =  strength  of  third  frequency  (47) 

xr5  =  double  strength  of  a  4-cycle  frequency 

xr6  =  xr4 

xr?  =  xr3 

xr8  =  xr2 

Thus,  a  Fourier  transform  of  N  points  can  handle  the  _N_  th 

2 

frequency  component  if  one  remembers  that  it  is  erroneously 
reported  at  double  strength.  If  one  desires  to  modify  the 
Fourier  components  of  a  signal  by  multiplying  them  by  image 
enhancement  weighted  factors,  the  two  elements  in  each  Pair, 
xr2  and  xrg,  xr3  and  xr7 ,  and  finally  xr4  and  xr6  must  be 


multiplied  by  the  same  factor.  If  the  upper  values  are 
ignored,  enhancement  is  incorrect. 

Remember,  before  the  application  of  the  inverse 
Fourier  transform,  the  results  must  be  divided  by  N  and 

2 

upon  return  to  physical  space  after  the  inverse  Fourier 
transform  is  applied,  the  results  must  be  divided  by  two 
for  proper  scaling. 

Subroutine  FT16H 

This  subroutine  invokes  subroutine  FT  to  do 
Fourier  transform  processing  of  the  image  data.  Each 
16-element  horizontal  picture  line  is  processed  in  turn 
from  the  top  to  the  bottom  of  the  16-by-16  input  image 
matrix. 

Subroutine  FT2D 

This  subroutine  takes  the  two-dimensional 
Fourier  transform  of  the  N-by-N  real  array  XR,  alters 
the  spatial  frequency  components  by  the  image  enhancement 
weighted  factors  that  the  user  furnishes,  and  then  takes 
the  inverse  Fourier  transform.  The  weighted  factors  are 
supplied  to  the  N  length  H  and  V  arrays  for  image  enhance- 
2 

ment  in  both  horizontal  and  vertical  directions.  In  this 
conventional  Fourier  transform  algorithm,  trigonometric 
values  are  computed  once  as  in  subroutine  FT  so  that  the 
number  of  required  trigonometric  calculations  is  reduced 
from  N4  to  N,  thus  saving  computer  time. 

Subroutine  FT16D2 

This  subroutine  collects  the  input  image  data 
and  applies  subroutine  FT2D  to  it  to  enhance  the  16-by-16 
image  array  via  the  two-dimensional  Fourier  transfoim 
algorithm . 

Fast  Fourier  Transform 

This  particular  algorithm  was  adapted  and 
debugged  for  Fortran  operation  from  an  algorithm  in  the 
"focal"  language.  It  was  originally  implemented  by  an 
investigator  at  the  Digital  Equipment  Corp.  in  Maynard, 
Massachusetts  (ref.  4).  Whereas  the  original  algorithm 
implemented  the  fast  Fourier  transform  only,  the  author 
has  modified  it  to  incorporate  image  enhancement  weighted 
factors  to  take  the  inverse  fast  Fourier  transform  and  to 


scale  the  data  upon  return  to  pixel  space.  The  user  is 
required  to  furnish  NU^ that  iSj the  log  to  the  base  2  of 
the  number  of  points  N.  Also,  the  array  XR  containing  N 
pixel  values  to  be  modified  and  the  array  E  containing  the 
N  weighted  factors  are  required.  Note  that  subroutine 
2 

REVERS  is  required  by  the  fast  Fourier  algorithm  to  properly 
compute  the  trigonometric  values. 

Fast  Fourier  Transform  Scans 

Subroutine  FFTH  and  FFTHV  invoke  the  fast  Fourier 
transform  algorithm  to  operate  on  the  input  target  image. 

The  first  scans  the  image  horizontally  from  top  to  bottom. 
The  second  scans  through  the  image  in  two  directions,  first 
across  horizontal  lines  from  top  to  bottom,  then  down 
vertical  columns  from  left  to  right  to  operate  on  the  data 
both  horizontally  and  vertically. 

Examples  of  Image  Enhancement  Demonstrations 

Figures  2,  3,  4,  and  5  are  examples  of  actual  image 
enhancement  incorporating  the  input  image  array  with 
background  value  19,  signal  value  20,  and  rms  noise  level  0. 

Figure  2  shows  digital  processing  in  the  horizontal 
direction,  one  complete  16-element  picture  line  at  a  time, 
by  the  fast  Fourier  transform  invoked  by  subroutine  FFTH. 

The  attempt  fails  catastrophically;  extensive  false  detail 
is  generated;  for  example,  the  gun  in  rows  3,  4,  and  5  of 
the  images  is  completely  destroyed.  One  can  speculate  that 
the  incapacity  of  the  fast  Fourier  transform  to  handle 
this  two-tone  image  is  due  to  the  limitations  of  the 
transform  bandwidth.  Where  N  =  16,  the  transform  can 
properly  sample  only  8  spatial  frequency  components,  a  dc 
term,  and  one  through  7  cycles  of  sinusoidal  oscillation 
for  each  data  field  length  of  16  points.  Therefore,  seven 
terms  are  not  enough  in  the  Fourier  series  expansions  of 
the  image  signal  to  accurately  represent  the  bandwidth 
needed  to  characterize  the  sharp  transitions  at  the  edges 
of  the  various  tank  target  elements. 

Figure  3  elici'.s  better  results  and  shows  digital 
processing  in  the  horizontal  direction  by  the  Haar  transform 
invoked  by  subroutine  HAR4H ,  i.e.,  four  elements  for  each 
transform  and  4  transforms  for  each  horizontal  picture 
line.  The  m-values,  0.5,  10  and  15,  are  the  enhancement 
parameters  and  weight  the  strengths  of  the  dc  term  and 
terms  of  one  and  two  cycles  of  square-wave  oscillation, 
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Figure  3.  Haar  transform,  4-element  horizontal  scan. 
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Figure  4.  Hadamard  transform,  horizontal 
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igure  5.  Haar  transform,  16-element  horizontal  scan. 


respectively,  for  each  data  field  length  of  four  pixels. 
Notice  that  contrast  is  improved  both  by  boosting  the 
signal  level  for  the  gun  and  the  other  tank  components 
through  the  second  two  of  the  three  m-values,  and  by 
reducing  the  background  almost  50  percent  through  the 
first  m- value  which  takes  out  some  of  the  dc  term. 

Figure  4  shows  horizontal  processing  by  the  Hadamard 
4-by-4  matrix,  through  subroutine  HAD4H ,  with  m-values  0.5, 
15,  10,  and  10.  The  first  governs  the  background;  the 
second,  the  highest  spatial  frequency;  the  third,  the 
next  highest,  and  so  on.  The  four  spatial  frequency 
components  are  not  as  systematically  arranged  as  in  the 
Haar  transform  because  the  sequencies  or  frequencies  of 
value  change  are  out  of  order  for  this  particular  algorithm 
for  the  Hadamard  transform. 

Finally,  figure  5  shows  processing  in  the  horizontal 
direction  by  subroutine  HAR16H,  i.e.,  via  the  HAAR  transform, 
taking  16  elements  or  each  line  of  pixels  at  a  time.  The 
M-values  are  0.5,  1,  5,  10,  and  15,  and  govern  the  d.c. 
term  and  terms  of  1,  2,  4,  and  8  cycles  of  spatial  frequency 
oscillation  for  each  data  field  length  of  16  pixel  points, 
respectively.  This  example  produces  the  best  contrast 
improvement  results  because  of  the  greater  flexibility 
afforded  by  the  larger  number  of  available  spatial  frequency 
enhancement  parameters  that  can  be  varied.  For  example, 
because,  in  figure  3,  the  Haar  4-element  scan  "sees"  the 
data  in  blocks  of  four,  it  never  sees  the  transition 
between  the  background  and  the  tread  of  the  tank.  Round 
off  in  the  integer  expression  of  the  output  obscures  the 
tread  by  blending  it  into  the  background,  giving  elements 
the  same  pixel  value  of  ten.  The  16-element  scan  of 
subroutine  HAR16H  does  see  the  tread  and  displays  it  with 
an  intensity  about  30  percent  higher  than  the  background. 

Note  again  that  in  these  scans,  where  contrast 
differences  are  pulled  out  or  enhanced  horizontally,  the 
results  over  the  scan  of  several  horizontal  lines  is  the 
enhancement  of  vertical  edges  or  lines  oriented  prefer¬ 
entially  in  the  vertical  direction. 

CONCLUSION 

An  interactive  program  with  the  capability  of  test  - 
demonstrating  the  effects  of  image  enhancement  and 
processing  algorithms  currently  exists  on  the  CDC  6500 
computer  system  at  the  ARRADCOM  Dover  site.  The  current 
status  of  the  program  exhibits  nine  horizontal,  vertical, 
and  two-dimensional  scanning  algorithms  employing  three 


transform  processing  techniques:  specifically  the  Haar, 
Fourier,  and  Hadamard  types  in  modular  form.  Other 
algorithms  of  interest  to  the  user  can  be  simply  tacked 
on  as  Fortran  callable  subroutines.  Then  algorithms  can 
employ  the  package's  capability  of  contrast,  and  signal- 
to-noise  measurement  of  input  imagery  data. 
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APPENDIX  A.  PROGRAM  IMAGE 


PROGRAM  IMAGE ( INPUT , OUTPUT , TAPE 4  =  6 5 , TAPE 61 =10 0 , TAPE 6 2=10 0 ) 
DIMENSION  IN  (16, 16) , IOUT (16,16) , IX ( 56 ) , IY  ( 56 ) 

DATA  IX/1641, 2986, 2986, 31 54, 31 54, 2658, 2662, 3094, 

&  31 38, 3146, 3162, 31 14, 3082, 30 3 8, 301 8, 3018, 24 90, 2490, 1981, 1981, 

&181 7, 181 7, 1469, 1469, 1641, 164 1,1 753, 1753, 1729, 2994, 2994, 3162, 
&3162, 2662, 266 2, 31 34, 3146, 3158, 31 30, 3082, 30 50, 30 14, 2998, 2518, 
&2494, 2494, 1977, 1977, 181 3, 1813, 1477, 1477, 1641, 1641, 1753, 1753/ 

DATA  I Y/l 991, 1991, 2071, 2071,2427,2427,2567,2811, 

&281 5, 2831, 287 5, 2883, 28 79, 2875, 2855, 2827, 2599, 2703, 2703, 2611, 
&2611, 2435, 2435, 2079, 20 79, 1991, 1991, 1991, 408, 404, 492, 492, 

&844, 848, 996, 1204, 1232,1268,1284,1288,1284,1272,1252,1020, 

&1020, 1116, 1116, 10 28, 1028, 84 8, 848, 488, 488, 408, 408, 40 8/ 

CALL  I N ITT ( 3 0 ) 

PRINT  * , " DO  YOU  WANT  TO  TEST  NUMBER  SEQUENSES  ON  YOUR  ALGORITHM?" 
CALL  BEEP(1,. 1,0.0) 

READ  50, NS 
FORMAT  (A3) 

IF (NS  . EQ.  "YES")  GO  TO  2 

PRINT  *, "ENTER  BACKGROUND,  SIGNAL  AND  NOISE  LEVELS.  " 

CALL  BEEP (3,  .  1,  .  2) 

READ  * , B , S , RMS  $CONTRA  =  100.*(S-B)/S 

IF  (RMS.EQ.0.)  STN=100.  $IF  (RMS. NE. 0.0)  STN=S/RMS 

PRINT  *, "SPECIFY  WHICH  OF  THE  FOLLOWING  TRANSFORMS  BY  NUMBER:" 

PRINT  * ,  "  1-HAR  (16,16)  MATRIX,  HORIZONTAL  SCAN." 

PRINT  * , " 2 -HAR  (16,16)  MATRIX,  VERTICAL  SCAN." 

PRINT  * , " 3-HAR (16,16)  MATRIX,  HORIZONTAL  THEN  VERTICAL  SCAN." 
PRINT  * , " 4-HAR (4,4)  MATRIX,  HORIZONTAL  SCAN." 


PRINT  *, "1-HAR  (16, 16)  MATRIX,  HORIZONTAL  SCAN." 
PRINT  * , " 2-HAR  (16,16)  MATRIX,  VERTICAL  SCAN." 

PRINT  * , "3-HAR (16,16)  MATRIX,  HORIZONTAL  THEN  VERT 
PRINT  *, "4-HAR  (4, 4)  MATRIX,  HORIZONTAL  SCAN." 
PRINT  * , " 5-H ADAMARD (4,4)  MATRIX,  HORIZONTAL  SCAN." 
PRINT  *, "6-FOURIER  TRANSFORM,  HORIZONTAL  SCAN." 
PRINT  * , " 7-FOURIER  TRANSFORM,  2-DIMENSIONAL  SCAN." 
PRINT  * , " 8  FFT  HORIZ.  SCAN" 

PRINT  * , " 9-FFT  HORIZ.  THEN  VERT.  SCAN" 

CALL  BEEP ( 1 , .1, .1) 

13  READ  *,  INDEX 

IF (INDEX  .LE.  9)  GO  TO  14 

PRINT  *, "INCORRECT  INCEX,  TRY  AGAIN.  " 

GO  TO  13 

14  IF (NS  .EQ.  "YES")  GO  TO  5 
DO  6  1=1,16 

DO  6J=1 , 16 

6  IN  (I , J) =3 
DO  7  1=1,3 

IN (2+1 , 14-1) =S 

7  IN ( 5 , 6+1 ) =S 
DO  8  1=6,7 
DO  8  J=6 , 10 

8  IN  (I , J ) =S 
DO  9  1=8,11 
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9  J=4, 13 
9  IN ( I ,  J )  =  S 
DO  10  J  =5 , 1 2 

10  IN ( 12 , J ) =S 

SQ  =  RMS  *SQRT ( (2.*16.+1.)/(6.*16.) ) 

IF (RMS  . EQ.  0.0)  C  =  0.0 
IF ( RMS  .NE.  0.0)  C  =  RMS *RMS/SQ* 2 . 0 
DO  11  1=1,16 
DO  11  J=1 , 16 

11  IN (I , J)=IN (I , J)  +  (RANF (9) -0.  5) *C  +  0. 5 

5  IF  (INDEX. EQ.l)  CALL  HAR16H ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 2)  CALL  HAR16 V ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 3)  CALL  HAR16H V ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 4)  CALL  H AR4H ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 5)  CALL  HAD4H ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 6)  CALL  FTl 6H ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 7)  CALL  FT16D 2 ( IN , IOUT , NS ) 

IP  (INDEX. EQ. 8)  CALL  FFTH ( IN , IOUT , NS ) 

IF  (INDEX. EQ. 9)  CALL  FFTHV ( IN , IOUT , NS ) 

41  I F ( NS  .EQ.  "YES")  GO  TO  42 
PRINT  43 ,8, S, RMS 

43  FORMAT  (1H  " bACKGR=  "F4.1/1H  "SIGNAL^  "F4.1/1H  ”RMS=  "F4.1) 
PRINT  44, CONTRA, STN 

44  FORMAT ( 1H  , "CONTRAST  IS : " , F 4 . 1 " % "/1H  ,”S/N  RATIO  IS:'\F4.1) 
CALL  APS IG  (TOUT) 

CALL  PRMS(IOUT) 

CALL  CTRAST  ( IOUT ) 

51  CALL  HOME 
CALL  ANMODE 
DO  52  1=1,16 

52  PRINT53,  (IN  (I  ,J)  ,J  =  1, 16) 

53  FORMAT ( 1 H  ,17X, 1613) 

PRINT  " 

DO  54  1=1,16 

54  PRINT53,  (IOUT  (I  ,J)  ,J  =  1, 16) 

CALL  M0VA3S (IX  (1)  ,  IY  (1) ) 

DO  55  I  =  2, 28 

5  5  CALL  DRWABS (IX  (I )  ,  IY  (I  )  ) 

CALL  MOVA3S (IX  (29)  , IY  (29) ) 

DO  56  I  =  30,56 
56  CALL  DRWABS (IX (I )  ,  IY  (I ) ) 

CALL  ANMODE 
CALL  BEEP ( 2 , . 1 , . 2) 

READ  * , DUMMY 

42  GO  TO  1 
END 

SUBROUTINE  8EEP(N,  T3 ,  TD) 

BELL  =  020702155555555555553 
3 LANK  =  555555555555555555558 
I  =  IFIX(10. 1 *TB ) 

J  =  IF IX  (5. 1 *TD ) 
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CALL  CONNEC (4,2) 

PRINT  " 

DO  1  K  =  1,N 
DO  2  L  =  1,1 

2  WRITE  (4)  BELL 

IF  (J.EQ.O)  GO  TO  1 
DO  3  M  =  1,J 

3  PRINT  4 , BLANK 
1  CONTINUE 

CALL  DI SCON (4,2) 

4  FORMAT  ( 1 H  + ,  A 1 ) 

RETURN 

END 

SUBROUTINE  APSIG(IOUT) 

DIMENSION  I OUT ( 16 , 16 ) 

PSIG=0 . 0 
DO  1  1=1,3 

PSIG=PSIG+IOUT (2+1 ,14-1) 

1  PSIG=P3IG+IOUT (5, 6+1 ) 

DO  2  1=6,7 

DO  2  J=6 , 10 

2  PSIG=PSIG+IOUT (I , J) 

DO  3  1=8,11 

DO  3  J=4 , 1 3 

3  PSIG=PSIG+IOUT (I , J) 

DO  4  J=5 , 12 

4  PSIG=PSIG+IOUT (12, J) 

PSIG=PSIG/64. 0 
PRINT5 , PSIG 

5  FORMAT ( 1H  ,"APSIG  IS  =  ",F6.2) 

RETURN 

END 

SUBROUTINE  PRMS(IOUT) 

DIMENSION  IOUT ( 16 , 16 ) 

CRMS=0 . 0 
DO  1  1=1,16 
DO  1  J=1 , 16 

1  CRMS=CRMS+IOUT (I , J) *IOUT (I , J) 

CRMS=SQRT (CRMS/25o. 0) 

PRINT2, CRMS 

2  FORMAT (1H  ,"PRMS  IS=  ",F6.2) 

RETURN 

END 

SUBROUTINE  CTRAST(IOUT) 

DIMENSION  IOUT  (16, 16)  ,R(16,16),A(12),G(12) 

DATA  A/"GUNH  =  "  ,  "GUNV="  ,  "COUH  =  "  ,  ,,COUV="  ,  " TURH  =  " 
5 " HULH  =  " , " HULV  =  " , "WHLH= " , "WHLV=" , ”TNKH= " , "TNKV= 
DO  2  I  =  1,16 
DO  2  J=l, 16 
R (I , J ) =IOUT (I , J) 

2  I F ( R ( I , J ) . EQ . 0 . )  R (I , J) =1 . E-14 


"TURV 

/ 
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G (1 ) = ( ( 2 . *R (3 , 13 ) -R ( 3 , 12 ) -R ( 3 , 14 ) ) /R ( 3 , 1 3 ) + (2*R (4 , 12) -R (4 , 11) -R (4 , 
113)  )  /R(4,  12)  +  (2*R (5,11 )-R (5,10) -R(5, 12) ) /R (5,11) ) /6 . 0*100.0 
G  ( 2 )  =  ( ( 2 *R  ( 3 , 13 ) -R { 2 , 13 ) -R { 4 , 1 3 ) ) /R ( 3 , 1 3 )  +  ( 2 *R (4 , 1 2 ) -R ( 3 , 1 2 ) -R ( 5 , 1 
12))/R(4,12)+(2*R(5,11)-R(4,11)-R(6,11))/R(5,11))/6.0*100.0 
G(3)=((R(5,7)-R(5,6))/R(5,7)  +  (R(5,9)-R(5,10))/R(5,9))/2.0*100.0 
G(4)  =  ((R(5,7)-R(4,7))/R(5,7)  +  (R(5,8)-R(4,8))/R(5,8)  +  (R(5,9)-R(4,9) 

I) /R(5, 9) )/3. 0*100.0 

G(5)=((R(6,6)-R(6,5))/R(6,6)  +  (R(6,10)-R(6,11))/R(6,10)  +  (R(7,6)-R(7 
l,5))/R(7,6)+(R(7,10)-R(7,ll))/R(7,10))/4. 0*100.0 
G(6)=( (R( 6,6) -R (5,6) ) /R( 6,6)+ (R( 6,10) -R( 5,10) )/R( 6, 10) )/2. 0*100.0 
G(7)=((R(8,4)-R(8,3))/R(8,4)+(R(8,13)-R(8,14))/R(8,13)+(R(9,4)-R(9 
1,3))/R(9,4)+(R(9,13)-R(9,14) )/R (9,13)+ (R (10,4) -R (10,3) )/R (10,4)+ (R 
2(10,13)-R(10,14))/R(10,13)+(R(11,4)-R(11,3))/R(11,4)+(R(11,13)-R(1 
31, 14) ) /R (11, 13) ) /8. 0*100. 0 

G ( 8 ) = ( ( R ( 8 , 4 ) -R ( 7 , 4 ) ) /R ( 8 , 4 ) + (R ( 8 , 5 ) -R ( 7 , 5 ) ) /R ( 8 , 5 ) + (R (8 , 11 ) -R ( 7 , 1 

II) ) /R ( 8 , 11 ) + (R ( 8 , 12 ) -R ( 7 , 12 ) ) /R ( 8 , 12 ) + (R ( 8 , 13 ) -R (7, 13) ) /R(8, 13) + (R 
2 (11, 4) -R (12, 4) )/R (11,4)+ (R (11,13) -R ( 12, 13) ) /R(ll, 13) ) /7 . 0*1 00.0 

G  (9 ) =  (  (R ( 1 2 , 5 ) -R ( 1 2 , 4 ) ) /R ( 1 2 , 5 )  +  (R ( 1 2 , 12 ) -R (12 , 1 3 ) ) /R ( 1 2 , 12 ) ) *50 . 0 
G (10)= ( ( R (12 , 5) -R ( 1 3 , 5) ) /R (12,5)+(R(12,6)-R(13,6))/R(12,6)+(R(12,7 
1)-R(13,7))/R(12,7)+(R(12,8)-R(13,8))/R(12,8)+(R(12,9)-R(13,9))/R(1 
22,9)+(R(12,10)-R(13,10))/R(12,10)+(R(12,ll)-R(13,ll))/R(12,ll)+(R( 
312, 12) -R (13, 12)) /R (12, 12 ))/8. 0*100.0 
G  ( 11 )  =  (G  ( 1 )  +G  ( 3  )  +G  ( 5  )  +G  ( 7  )  +G  (9  )  )  /5  . 

G  ( 1 2  )  =  (G  ( 2 )  +G  ( 4  )  +G  (6  )  +G  (8  )  +G  (1 0  )  )  /5  . 

DO  3  1=1, 12 
3  PRINT  5 , A ( I )  , G  ( I ) 

5  FORMAT ( 1H  , A5 , F 7 . 2 , " %  "  ) 

RETURN 

END 

SUBROUTINE  HAAR (Q , E , NU ) 

DIMENSION  T ( 1 6 ) , Q ( 3 0 )  ,E(5) 

N  =2**NU 
NO  =N 
ME  =  NU+2 
IP  =  (-2 ) *N 
IR  =  0 

5  NO  =  NO/2 

IP  =  IP  +  4  *NQ 
IR  =  IR  +  2  *NO 
ME  =  ME -1 
DO  6  J  =1 ,  NO 
L  =  2  *J 
K  =  L-l 
NO J  =  NO+J 

T(NOJ)  =  E (ME) * (Q (IP  +  K) -Q (IP+L)  ) 

6  Q (IR+J)  =Q (IP+K) +Q(IP+L) 

IF  (NO.GT.2)  GO  TO  5 

T ( 2 )  =E (2) * (Q (IR+1 ) -Q (IR+2) ) 

T  ( 1 )  =E(l)*(Q(IR  +  l)+Q(IR  +  2)  ) 

3(1)  = (T (1 ) +T (2) ) /FLOAT  (N) 

3(2)  = (T (1 ) -T (2) ) /FLOAT  (N) 
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IP=  -1 
IR  =  0 

7  N  2  =  NO 

X  =  FLOAT (N 2) /FLOAT (N) 

NO  =  NO *2 

IP  =  IP  +  NO/4 

IR  =  IR  +  NO/2 

DO  8  J  =1 ,  NO 

ND  =  (J+l)/2 

ID  =  N2  +  ND 

Y  =  (-1.0)**(J+1)*X 

8  Q  ( IR+J )  =  Q (IP+ND) +Y*T (ID) 

IF  (NO.LT.N)  GO  TO  7 
RETURN 

END 

SUBROUTINE  HAR16H ( IN ,  IOUT,  NS) 

DIMENSION  IN (16, 16) , IOUT (16, 16) ,Q(30) ,E(5) 
PRINT  *," ENTER  5  M-VALUES." 

CALL  BEEP ( 5 , . 1, . 1) 

READ  *,  (E  (J)  , J-1,5) 

CALL  ERASE 

PRINT  * , "HAR16H  5HOR.  ENH" 

PRINT  1,  (  ( J , E ( J ) )  , J  =  l, 5) 

1  FORMAT  ( 5  (1H  ,"M",I1,"  =  ",F4.1,/)) 

NU  =  4 

M  =  14 

IF  (NS  .EQ.  "YES")  GO  TO  2 
DO  3  I  =  1,16 
DO  4  J  =  1,16 

4  Q(J)  =  IN  (I,J) 

CALL  HAAR (Q , E , NU ) 

DO  5  K  =  1,16 

5  IOUT (I ,  K )  =  Q (M+K ) +0 . 5 
3  CONTINUE 

GO  TO  6 

2  PRINT  *, "ENTER  16  INPUTS." 

CALL  BEEP (16 , . 1, .2) 

READ  *, (Q ( J ) , J=1 , 16 ) 

CALL  HAAR (Q , E , NU ) 

PRINT  *,  ( Q ( M  + J )  , J=1 , 16 ) 

6  RETURN 
END 

SUBROUTINE  HAR16V(IN,  IOUT,  NS) 

DIMENSION  IN (16, 16) ,IOUT(16,16) ,Q(30) ,E(5) 
IF (NS  .EQ.  "YES")  GO  TO  2 
PRINT  *, "ENTER  5  M-VALUES." 

CALL  BEEP ( 5 , . 1, . 1) 

READ  *,  (E  (J)  ,  J  =  1,5) 

CALL  ERASE 

PRINT  * ,  "HAR16V  5VER.  ENH" 

PRINT  1,  ( ( J , E  ( J ) )  f J-1,5) 
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1  FORMAT  ( 5 ( 1H  , " M " , 1 1 ,  "  =  " ,F4.1/)) 

MLJ  =  4 

M  =  14 

DO  3  I  =  1,16 
DO  4  J  =  1,16 

4  Q(J)  =  IN  ( J ,  I ) 

CALL  H AAR (Q , E , NU ) 

DO  5  K  =  1,16 

5  IOUT {K,I)=Q(M+K)+0. 5 

3  CONTINUE 

2  RETURN 
END 

SUBROUTINE  HAR16HV ( IN , IOUT , NS ) 

DIMENSION  IN  (16, 16)  , IOUT  (16, 16)  ,Q(30),H(5),V(5) 

IF (NS  . EQ.  "YES”)  GO  TO  1 

PRINT  *, "ENTER  5  HORIZONTAL  M-VALUES." 

CALL  BEEP  (5, . 1, . 1) 

READ  *,  (H  (J)  , J  =  l, 5) 

PRINT  *, "ENTER  5  VERTICAL  M-VALUES" 

CALL  BEEP ( 5 , . 1, . 1) 

READ  *, (V(J) ,J=1, 5) 

CALL  ERASE 

PRINT  * , " HAR 16HV  5HOR.  ENH" 

PRINT  2, ( ( J , H ( J ) ) , J=1 , 5) 

PRINT  * , " 5  VERT  ENHANCERS" 

PRINT  2, ( (J,V(J) ) , J*l, 5) 

2  FORMAT  ( 5  ( 1H  ,"M",I1,”  *  ”,F4.1,/)) 

M=14 

MU  =  4 

DO  3  I  =  1,16 

DO  4  J  =  1,16 

4  Q ( J )  =  IN  (I,J) 

CALL  HAAR (Q , H , NU ) 

DO  5  K  =  1,16 

5  IOUT (I ,K)  =  Q (M+K) +. 5 

3  CONTINUE 

DO  6  I  =  1,16 

DO  7  J  =  1,16 

7  Q  ( J )  =  IOUT  (J, I) 

CALL  HAAR (Q , V, NU ) 

DO  8  K  =  1,16 

8  IOUT  (K, I)  =  Q  (M+K) +. 5 

6  CONTINUE 
1  RETURN 

END 

SUBROUTINE  HAR4H ( IN ,  IOUT,  NS) 

DIMENSION  IN (16,  16) , IOUT (16, 16),Q(6),E(3) 

IF  (NS. EQ. "YES”)  GO  TO  6 
PRINT  *,  "ENTER  3  M-VALUES." 

CALL  BEEP ( 3 , . 1 , . 1 ) 

READ  *,  (E  (J)  , J  =  l, 3) 
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CALL  ERASE 

PRINT  * , " HAR4H  5H0R.  ENH" 

PRINT  1, ( (J,E (J) ) ,J=1, 3) 

I  FORMAT  ( 3 ( 1H  ,"M",I1,"  =  ",F4.1,/)) 

NU  =  2 
M  =  2 

DO  3  I  =  1,16 

DO  3  J  =  1,4 

DO  4  K  =  1,4 

4  Q(K)  =  IN  (1 , 4*  ( J-l )  +K) 

CALL  HAAR{Q,E,NU) 

DO  5  L  =  1,4 
K  =  4  * ( J-l ) +L 

5  IOUT (I ,K)  =  Q (M+L) +. 5 
3  CONTINUE 

6  RETURN 
END 

SUBROUTINE  HAD (Q,E,NU) 

DIMENSION  T (4) ,Q(4) ,E (4) 

N  =  2**NU 
IT  =  0 

1  IT  =  IT  +1 

T  ( 1 )  =  Q  ( 1 )  +  Q  (2)  +  Q  (3  )  +  Q  (4  ) 

T  ( 2 )  =  Q  ( 1 )  -  Q  ( 2 )  +  Q  (3  )  -  Q(4) 

T(3)  =  Q(l)  +  Q(2)  -  Q(3)  -  Q(4) 

T  (4  )  =  Q  ( 1 )  -  Q  ( 2)  -  Q  (3  )  +  Q(4) 

DO  2  I  =  1,N 

IF (IT  .EQ.  1)  T ( I )  =  T ( I )  *  E ( I ) 

IF (IT  . EQ.  2)  T (I)  =  T  (I)  /  N 

2  Q  (I )  =  T  ( I ) 

IF (IT  .LT.  2)  GO  TO  1 

RETURN 

END 

SUBROUTINE  HAD4H (IN ,  IOUT,  NS) 

DIMENSION  IN (16, 16) , IOUT (16, 16) ,Q (4) ,E(4) 
NU  =  2 
N  =  2**NU 
M  =  0 

PRINT  *, "ENTER  ",N,"  HORIZONTAL  ENHANCERS. 
CALL  BEEP (N , . 1, .1) 

READ  *, (E (J) , J=1,N) 

CALL  ERASE 

PRINT  * , "HAD 4 6  EHOR.  ENH" 

PRINT  1, ( ( J , E ( J ) ) ,J=1,N) 

1  FORMAT  ( 4  ( 1H  ,"M",Il,"  =  ",F4.1,/)) 

IF (NS . EQ. "NO")  GO  TO  2 
DO  3  I  =  1,16 

DO  3  J  -  1, N 

DO  4  K  =  1,N 

4  Q  (K )  =  IN  ( I ,  N*  ( J-l )  +  K) 

CALL  HAD (Q , E , NU ) 
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DO  5  L  =  1,N 
K  =  4* ( J-l )  +  L 

5  IOUT (I  ,K)  =  Q (M-t-L) +  . 5 
3  CONTINUE 

GO  TO  6 

2  PRINT  *, "ENTER  " ,N,"  INPUTS.  " 

CALL  BEEP (N  ,  .  1 ,  .  1 ) 

READ  *,  (Q ( J )  , J  =  1 , N) 

CALL  HAD (Q ,  E , NU ) 

PRINT  *,  (Q  (M+J)  , J=1,N) 

6  RETURN 
END 

SUBROUTINE  FT(XR,NU,E) 

DIMENSION  CO ( 16 ) ,SI (16) , XR (16) ,XI (16) ,  E (8) ,FR(16) ,FI (16) 
N  =  2**NU 
NP  =  N/2 

TPN  =  8.  *  ATAN ( 1 . )  /  N 
DO  1  I  =  1,N 
X  =  TPN  *  (1-1) 

CO (I)  =  COS  (X) 

1  SI (I)  =  SIN  (X) 

IC  =  0 

6  IC  =  IC  +  1 

SIGNR  =  FLOAT ( (-1) **IC) 

SIGNI  =  -1.  *  SIGNR 

IF ( IC  .EQ.  1)  SCALE  =  NP 

IF ( IC  .EQ.  2)  SCALE  =2.0 

DO  2  I  =  1,N 

FR ( I }  =  FI  (I)  =0. 

DO  2  J  -  1,N 

Ml  =  ( I -1 )  *  (J-l) 

M2  =  MOD (Ml , N)  +  1 

FR  (I )  =  FR  ( I )  +  CO (M2 ) *XR ( J )  +  3IGNR*SI (M2) *XI (J) 

2  FI ( I )  =  F I ( I )  +  CO (M2) *XI (J)  +  SIGNI*SI (M2) *XR(J) 

DO  3  I  *  1,N 

XR ( I )  =  FR ( I )  /  SCALE 

3  XI (I)  =  F I ( I )  /  SCALE 
IF ( IC  .EQ.  2)  GO  TO  4 
DO  5  I  *  1 ,  NP 

J  =  N+2  -I 

XR ( I )  =  XR(I)  *  E  (I) 

XI (I)  =  XI  (I)  *  E  (I) 

IF ( I  .GT.  1)  XR ( J )  =  XR ( J )  *  E  (I ) 

5  I F ( I  .GT.  1)  XI(J)  =  XI  (J)  *  E  ( I ) 

GO  TO  6 

4  RETURN 
END 

SUBROUTINE  FT16H ( IN , IOUT, NS ) 

DIMENSION  IN (16, 16) , IOUT (16, 16) ,XR(16) ,E(8) 

NU  =  4 
N  =  2**NU 
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NP  =  N/2 

IF (NS  . EQ.  "YES")  GO  TO  1 

PRINT  *, "ENTER  ",NP,"  HORIZONTAL  ENHANCERS.  " 

CALL  3EEP (NP, . 1, . 1) 

READ  *,  (E  (J) , J=1,NP) 

CALL  ERASE 
DO  2  1=1, N 
DO  3  J=1 , N 

3  XR (J)=IN ( I ,  J ) 

CALL  FT (XR, NU ,  E) 

DO  4  K=1,N 

4  IOUT (I , K) =XR (K) +. 5 

2  CONTINUE 

PRINT  * , "FT.  " , NP , "  HOR.  ENH" 

PRINT  5,  ( ( J , E  ( J ) )  , J  =  1 , NP ) 

5  FORMAT  (8 (1H  ,"M",Ilf"  =  ",F4.1,/)) 

1  RETURN 

END 

SUBROUTINE  FT2D (XR, NU , H , V) 

DIMENSION  CO (16) ,SI (16) ,XR(16,16) , XI (16,16) ,H(8) ,V(8) , 
&FR (16,16) ,  FI (16,16) 

N  =  2**NU 
NP  =  N/2 

TPN  =  8.  *  ATAN ( 1 . )  /  N 
DO  1  1  -  1, N 
X  =  T  PN  *  ( I  - 1 ) 

CO(I)  =  COS(X) 

1  31  (I)  =  SIN  (X ) 

IC  =  0 

3  IC  =  IC  +  1 

SIGNR  =  FLOAT ( (-1 ) **IC) 

SIGNI  =  -1.0  *  SIGNR 

IF ( IC . EQ.  1)  SCALE  =  N*N/4 

IF ( IC  .EQ.  2)  SCALE  =  4.0 

DO  4  L  =  1,N 

DO  4  K  =  1,N 

FR(L,K)  =  FI  (L, K)  =  0. 

DO  4  J  -  1, N 
DO  4  I  -  1, N 

Ml  =  (L-l )  *  (J-l)  +  (K-l)  *  (1-1) 

M2  =  MOD (Ml , N)  +1 

FR (L , K )  =  FR (L,K) +CO(M2) *XR(J , I) +SIGNR*SI (M2) *XI ( J, I ) 

4  FI (L,K) =FI (L , K) +CO (M2 ) *XI ( J , I ) +SIGNI *SI (M2) *XR(J,I) 

DO  5  I  -  If N 

DO  5  J  =  1 , N 

XR (I , J)  =  FR ( I , J )  /  SCALE 

5  XI  (I,J)  =  FI  (I , J)  /  SCALE 
IF  ( IC  .EQ.  2)  GO  TO  6 

DO  7  I  -  1,  NP 
K  =  N+2  -I 
DO  7  J  -  If N 
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XR(J,I)  =  XR ( J , I )  *  H { I ) 

XI  (J,I)  =  XI (J,I)  *  H (I) 

XR ( I , J )  =  XR ( I , J )  *  V(I) 

XI  (I,J)  =  XI (I,J)  *  V(I) 

IF(I.GT.l)  XR ( J , K)  =  XR ( J, K)  *  H(I) 

IF(I.GT.l)  XI ( J , K)  =  XI { J ,  K )  *  H ( I ) 

I F  ( I .  GT .  1)  XR  (K ,  J )  =  XR  (K ,  J)  *  V(I) 

7  IF(I.GT.l)  XI (K ,  J )  =  XI ( K , J )  *  V(I) 

GO  TO  3 
6  RETURN 
END 

SUBROUTINE  FT16D2(IN,  IOUT,  NS) 

DIMENSION  IN (16, 16) , IOUT (16, 16) , XR (16,16) ,H(8) ,V(8) 

NU  =  4 
N  =  2**NU 
NP  =  N/2 

IF (NS. EQ. "YES")  GO  TO  1 

PRINT  *, "ENTER  " , NP , "  HORIZONTAL  SPATIAL  FREQUENCY  ENHANCERS.  " 
CALL  BEEP (N ,  .1,  .1) 

READ  *,  (H  (J)  ,  J  =  1 , NP ) 

PRINT  *, "ENTER  " ,NP,"  VERTICAL  SPATIAL  FREQUENCY  ENHANCERS.  " 
CALL  BEEP (NP,  .1,  .1) 

READ  *, (V(J) , J=1,NP) 

DO  2  I  =  1,N 

DO  2  J  =  1,N 

2  XR ( I , J )  =  IN (I, J) 

CALL  FT2D(XR,NU,H,V) 

DO  3  I  =  1,N 

DO  3  J  =  1,N 

3  IOUT ( I , J )  =  XR ( I , J ) + . 5 
CALL  ERASE 

PRINT  * , "FT 2D  ",NP,"  HOR.  ENH" 

PRINT  10, ( (J,H(J) ) ,J=1,NP) 

10  FORMAT ( 8 ( 1H  ,"M",Il,"  =  ",F4.1,/)) 

PRINT  * , "FT 2D  ",NP,"  VERT.  ENH" 

PRINT  10, ( ( J  ,  V  ( J ) ) ,J=1,NP) 

1  RETURN 
END 

SUBROUTINE  FFT (XR, NU , E ) 

DIMENSION  XR(16)  ,XI (16) ,E  (8) 

INTEGER  S,Q,H,P,C,U 
IC=  0 
N  =  2**NU 

TPN  =  8.  *  ATAN ( 1 . )  /  FLOAT (N) 

X=3. O/FLOAT (N) 

DO  23  1=1, N 

23  XI (I) =0. 0 

24  IC=IC+1 

SIGNR  =  FLOAT ( (-1 ) **IC) 

SIGNI  =  -1.  *  SIGNR 
H=1-NU 
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tr 


I 


{ 


L=1 

NP=N/2 

Q=NP 

S=NP 

2  XSTOR=XR(Q) +XR(Q+S) 

XR  (Q+S )  =XR  (Q)  -XR  (Q+S ) 

XR(Q)=XSTOR 

IF  (IC.EQ. 1)  GO  TO  3 
XSTOR=XI (Q) +XI (Q+S) 

XI (Q+S) -XI (Q ) -XI (Q+S) 

XI (Q) =XSTOR 

3  IF  (Q.LE. 1)  GO  TO  5 

4  Q=Q-1 
GO  TO  2 

5  IF  (L.EQ.NU)  GO  TO  14 

6  L=L+1 
S=S/2 
H=H+1 
P=N 

7  C=1 

8  U=FLOAT (P-1 ) *2. 0**FLOAT (H) +X 
CALL  REVERS (U ,  K,  IR,  NU ) 

CO-COS (TPN*K) 

SI  =  SIN  (TPN*K ) 

GR  =  CO  *  XR (P)  +SIGNR  *  SI  *  XI (P) 
GI  =  CO  *  XI (P)  +  SIGNI  *  SI  *  XR (P ) 
Q=P-S 

XR ( P ) =XR (Q ) +GR 
XR  (Q  )  =XR  (Q)  -GR 
XI (P)=XI (Q) +GI 
XI (Q ) =XI (Q)-GI 

9  P=P-1 

10  IF (C. EQ. S)GO  TO  12 

11  C=C+1 

GO  TO  8 

12  IF (P . EQ. S ) GO  TO  5 

13  P=P-S 
GO  TO  7 

14  Q=N 

15  P=P-1 
Q=Q-1 

CALL  REVERS (Q,P, IR,NU) 

Q=Q  +  1 
P=P  +  1 

16  IF  (P.GE.Q)  GO  TO  18 

17  XSTOR=XR (P) 

XR (P) =XR (Q) 

XR (Q) =XSTOR 
XSTOR-XI (P) 

XI (P ) =XI (Q) 

XI (Q ) -XSTOR 
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18  IF(Q.EQ.l)  GO  TO  27 

19  Q=Q-1 

GO  TO  15 

27  IF  (IC.EQ.l)  SCALE  =  NP 
IF  (IC.EQ.2)  SCALE  =2 . 0 
DO  28  I  =1 ,  N 
XR(I)=XR(I) /SCALE 

28  XI (I)=XI (I)/SCALE 

IF  (IC.EQ. 2)  GO  TO  29 
DO  30  1=1, NP 
J  =  N+2  -I 
XR ( I )  =  XR ( I )  *  E  (I) 

XI  (I)  =  XI  (I)  *  E  (I) 

IF(I.GT.l)  XR ( J )  =  XR { J )  *  E  ( I ) 

30  IF(I.GT.l)  XI(J)  =  XI(J)  *  E  ( I ) 

GO  TO  24 

29  RETURN 
END 

SUBROUTINE  REVERS ( IU , K , IR , NU ) 

DIMENSION  IR  (NU ) 

J  =  IU 
K  =  ND  =  0 
DO  1  1=1, NU 

1  IR ( I ) =0 

2  ND=ND+1 
IQ=J/2 

IR (ND) =J-2  *IQ 
J  =  IQ 

IF  (J.GT.O)  GO  TO  2 
DO  3  1 =i , NU 

3  K=K+IR (I ) *2** (NU-I ) 

RETURN 

END 

SUBROUTINE  FFTH ( IN , IOUT , NS ) 

DIMENSION  IN (16, 16) , IOUT (16,16) ,XR(16) ,E(8) 
NU  =  4 
N  =  2**NU 
NP  =  N/2 

I F (NS  .EQ.  "YES")  GO  TO  1 

PRINT  * , "ENTER  ",NP,"  ENHANC.  PRAMS." 

CALL  BEEP (NP,.l,.l) 

READ  *, (E (J) , J=1,NP) 

CALL  ERASE 
DO  2  1=1, N 
DO  3  J=1,N 

3  XR ( J ) =IN (I , J) 

CALL  FFT (XR, NU , E ) 

DO  4  K=1 , N 

4  IOUT ( I , K) =XR (K ) +.  5 
2  CONTINUE 

PRINT  * , " FFT  " ,NP,"  HOR.  ENH" 
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PRINT  6, ( ( J , E  (J) )  ,  J  =  1 , NP) 

6  FORMAT ( 8 ( 1H  ,"M",Il,"  =  ",F4.1,/)) 

1  RETURN 
END 

SUBROUTINE  FFTHV ( IN ,  IOUT,  NS) 

DIMENSION  IN (16, 16) , IOUT (16, 16) ,XR(16),H(8),V(8) 
NU  =  4 
N  =  2**NU 
NP  =  N/2 

IF (NS  . EQ.  "YES")  GO  TO  1 

PRINT  *, "ENTER  ",NP,"  HORIZONTAL  ENHANCERS.  " 
CALL  BEEP (NP,  0.1,  0.1) 

READ  *, (H (J) , J=1,NP) 

PRINT  * , "ENTER  " , NP , "  VERTICAL  ENHANCERS.  " 

CALL  BEEP  (NP ,  0.1,  0.1) 

READ  *, (V (J) , J=1,NP) 

CALL  ERASE 
DO  2  I  =  1 ,  N 
DO  3  J  =  1,N 

3  XR ( J ) =IN ( I , J ) 

CALL  FFT (XR, NU , H ) 

DO  4  K  «  1,8 

4  IOUT (I ,K) =XR (K) +. 5 

2  CONTINUE 

DO  5  I  »  1,8 
DO  6  J  =  1 , N 

6  XR ( J ) =I0UT ( J , I ) 

CALL  FFT (XR, NU , V) 

DO  7  K  =1 , N 

7  IOUT  (K ,  I )  =  XR  (K )  +.  5 

5  CONTINUE 

PRINT  * , "FFT  HORIZ." 

PRINT  9,  ( ( J , H  { J ) )  , J=1 , NP ) 

9  FORMAT  (8 (1H  ,"M",I1,"  =  ",F4.1,/)) 

PRINT  * , "FFT  VERT." 

PRINT  9, ( (J,V(J)) , J=1 , NP ) 

1  return 

END 


APPENDIX  B.  PROGRAM  FIG1 


PROGRAM  FIG1 ( INPUT , OUTPUT , TAPE 4  =  6 5 , TAPE 61  =  100 , TAPE 62  =  1 00 ) 
DIMENSION  IX(56) ,IY (56) ,  IN  (16,16) 

DATA  IX/1641, 2986, 2986, 3154,3154, 2658, 2662, 3094, 

&  31 38, 3146, 3162, 3114, 30  82, 303  8, 3018, 3018, 2490, 2490, 1981, 1981, 
&  181 7, 1817, 1469, 1469, 1641, 1641, 17  53,1753, 172  9,  2994,2994,3162, 
&  3162, 266  2, 266  2, 3134,3146, 3158, 3130, 3082, 30  50, 3014, 2998, 2  518, 
&  2494, 2494 ,1977,1977,1813,1813,1477,1477,1641,1641,17  53,17  53/ 
DATA  IY/1 991, 1991, 2071, 20 71, 24 27, 2427, 2567, 2811, 

&281 5, 2831, 287  5, 2883,2879, 287  5, 2855, 2827, 2599, 2703, 270  3, 2611, 
&261 1,243 5, 24 35, 2079, 2079, 199 1,1 991, 1991, 408, 404, 492, 492, 
&844, 848, 996, 1204,1232, 1268, 1284, 1288, 1284, 1272,1252, 1020, 
&1020, 1116, 1116, 1028, 1028, 848, 848, 438, 488, 408, 408, 408/ 

CALL  INITT  (30 ) 

DO  6  1=1,16 
DO  6J  =  1 , 16 

6  IN  (I , J) ="B" 

DO  7  1=1,3 

IN  (2  +  1 , 14-1 )  =  "G" 

7  IN ( 5 , 6+1 ) = "C" 

DO  8  1=6,7 

DO  8  J=6 , 10 

8  IN (I , J) ="T" 

DO  9  1=8,11 
DO  9  J=4 , 13 

9  IN (I,J)="H" 

DO  10  J=5 , 12 

10  IN (12, J)="W" 

51  CALL  HOME 
CALL  ANMODE 
DO  52  1=1,16 

52  PRINT53, (IN (I , J) , J=l, 16) 

53  FORMAT ( 1H  ,17X,16A3) 

CALL  M0VA8S ( I X ( 1 ) , IY (1 ) ) 

DO  55  I  =  2,28 

55  CALL  DRWABS ( I X ( I ) , IY (I) ) 

CALL  MOVABS (IX(29),IY(29)) 

DO  56  I  =  30,56 

56  CALL  DRWABS (IX(I),IY(I)) 

CALL  ANMODE 

END 


47 


DISTRIBUTION  LIST 


Director 

US  Army  Research  Office 
ATTN:  Library 

P.O.  Box  12211 

Research  Triangle  Park,  NC  27709 
Director 

Department  of  Defense 
Advanced  Research  and  Projects 
Agency 

Washington,  DC  20301 
Commander 

US  Army  Materiel  Development  & 
Readiness  Command 
ATTN :  DRCSA-R 

DRCSA-PP 

Washington,  DC  20315 
Commander 

US  Army  Aviation  Systems  Command 
12th  &  Spruce  Streets 
St.  Louis,  MO  63166 

Commander 

US  Army  Tank-Automotive  Research  & 
Development  Command 
ATTN:  DRDTA ,  Technical  Library 

Warren,  MI  48090 

Commander 

US  Army  Electronics  Command  & 
Devices  Laboratory 
Ft.  Monmouth,  NJ  07703 

Commander 

US  Army  Electronics  Command 
Atmospheric  Science  Laboratory 
White  Sands  Missile  Range 
ATTN:  DRSEL-BL-MS 

Mr.  Richard  Gomez 
White  Sands,  NM  88002 


49 


Commander 

US  Army  TECOM 

ATTN:  AMSTE-TA-A  (2? 

STEAP-TL 

STEAP-DS-LP 

Aberdeen  Proving  Ground,  MD  21005 
Commander 

US  Army  Harry  Diamond  Laboratories 
ATTN:  AMXDO-TD/002 

AMXDO-TIB 

2800  Powder  Mill  Road 
Adelphi,  MD  20783 

Commander 
US  Army  ARRADCOM 

ATTN:  DRDAR-SC  COL  A.  J.  Larkins 

DRDAR-SC  Dr.  D.  A.  Gyorog 
DRDAR-SCF ,  L.  Berman 
DRDAR-SCF-I ,  J.  Lehman 
DRDAR-SCF-IO ,  J.  Lester 
DRDAR-SCF- 10 ,  G.  Sivak  (12) 
DRDAR-SCP,  J.  Glennon 
DRDAR-MSA,  R.  Isakower 
DRDAR-TSS  (5) 

Dover,  NJ  07801 

r- c c  ^  ,T,firbr>  -j  I  Tn  forma on  Center 

Cameron  Scat  ion 
Alexandria,  VA  22314 

Commander 

US  Army  Armament  Materiel  and 
Readiness  Command 
ATTN:  DRDAR-LEP-L 

Rock  Island,  IL  61299 


o  i  l  C  w  w  i. 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SL  (Tech  Lib) 

White  Sands  Missile  Range,  NM  88002 

US  Armv  Materiel  Systems 
Analysis  Activity 

ATTN:  DRXSl’-MP 

;  b,-,rdor>n  Pre  v  i  ra  Orourr’ 


*  rr> 


1  n  c; 


Weapon  System  Concept  Team/CSL 
ATTN :  DRDAR-ACW 

Aberdeen  Proving  Ground,  MD  21010 

Technical  Library 

ATTN:  DRDAR-CLJ-L 

Aberdeen  Proving  Ground,  MD  21010 

Technical  Library 

ATTN:  DRDAR-TSB-S 

Aberdeen  Proving  Ground,  MD  21005 

Technical  Library 
ATTN:  DRDAR-LCB-TL 

Benet  Weapons  Laboratory 
Watervliet,  NY  12189 

Dr.  Thomas  S.  Huang 
Institute  for  Image  Technology 
1000  North  Western  Ave . 

West  Lafayette,  IN  47906 


